Exploring Rasters

CPLN-503 Modeling Geographic Objects

Our lab today will go through some common operations, and we will end with an estimate of the number of cells that meet the following set of criteria:

• Be within two miles of a road on which you can travel more than 50 miles per hour

• Be at an elevation between 300 and 600 feet

• Be within forest land

• Be a park

This kind of workflow is a fairly common on with raster analysis and can be thought of as a kind of a suitability analysis. I won’t link to any here, but google about and you will find many examples and tools of weighted overlay/suitability analyses using raster data. And, after today, you all will know enough to just be dangerous.

Behold! GIS PUSHEEN!

It’s Raster Time

We will be making use of a new package today: terra. terra is the successor package to the venerable raster package. raster is still operable (I think), but is slated for deprecation so if you are searching for further info or tutorials please make sure to double check what package is being used. After loading our packages we will read in the layers we need for today.

if(!require(pacman)){install.packages("pacman"); library(pacman)}
p_load(terra, sf)

elev <- rast(x = here::here("labs_exercises/Nov9_OR 2023-10-25_RasterIntro_DRAFT/Nov_9th In Class/Nov_9th In Class/Nov_9th/Raster_Intro/elevation/"))

roads <- st_read(here::here("labs_exercises/Nov9_OR 2023-10-25_RasterIntro_DRAFT/Nov_9th In Class/Nov_9th In Class/Nov_9th/Raster_Intro/Lancaster_Roads.shp"), quiet = TRUE)

Okay, we should see something new here. First, we are using the rast() function from terra to read in our raster file. This should feel familiar, BUT look at the file path here. We are pointing rast() to a folder, and not a file. This is because while rasters are simply image grids, they still require all of the necessary metadata and geographic information to be useful for our analyses. All of that information is stored in the directory that holds our images, and rast references that data so that we can actually make use of our rasters.

Grabbing Our Roads

Next, let’s choose all roads with speeds that are greater than or equal to 50mph. Considering that our roads dataframe is a sf object, what are some ways you could do this operation (note if you want to use dplyr remember to load in that package)?

roads_50 <- roads[roads$AVG_SPEED >= 50, ]

roads_50 <- vect(roads_50)

Before you view my approach, have you tried one on your own? Let’s break this down. I use a base-R approach to filter my roads dataframe for all those roads with an average speed greater than or equal to 50mph. Next, I take this filtered sf dataframe and place into a vect() call. vect() is a function from terra to read vector data files. We could have used it before to read in our roads shapefile, but if you are more comfortable using sf and then converting that is also an option.

Measuring Distances

Next, we want to get distance calculations from our roads in relation to our elevation raster. In other words, we want to know the distance to the nearest road for every pixel from our elevation raster and then we want to find what cells are within 2 miles of a road.

roads_dist <- distance(x = elev, y = roads_50)

|---------|---------|---------|---------|
=========================================
                                          

roads_2mi <- roads_dist <= 3218.69

Now, let’s take a look at our roads_dist raster. Notice that each cell has a unique value for the distance measure it has to the nearest roads feature.

Ideally, this should make sense to us. We run the distance() function with elev as our first argument and our roads_50() raster. Next we want to identify which cells in our distance raster are less than or equal 2 miles. I dn not quite understand why because our CRS is in feet, but distance() default behavior is to measure linear units in meters. I need to double check what’s happening there, but 3218.69 is a conversion for 2 miles to meters.

Working with Elevation

Next, we want to isolate those cells that are less than or equal to 300 feet in elevation and those that are greater than or equal to 600 feet in elevation.

elev_sub <- elev >= 300 & elev <= 600

plot(elev_sub)

Getting our new elev_sub layer is relatively straightforward as we can apply conditional statements to a raster directly using map algebra. The next step here will be a little more complicated, but bear with me. While we have our binary raster we also have a lot of NA cells. We’d like to convert those NA cells to 0s, and we also want to convert our FALSE cells to 0 and our TRUE cells to 1s. In order to do so we’ll make use of the classify() function and we will make a classification matrix.

If that sounds confusing, don’t worry. We’ll go through it step-by step.

elev_mat <- c(NA, NA, 0,
              1, 1, 1,
              0, 0, 0)

elev_mat <- matrix(data = elev_mat, byrow = TRUE, ncol = 3)

First, we will set up the structure of our matrix by creating a new vector called elev_mat. Then we feed elev_mat into the matrix() function, tell it to create our matrix by organizing our vector by row and to place it in 3 columns. So, why do this? In order to make use of the classify() function it needs to know what values to classify and what to classify them into. In order to do that we create a matrix with three columns. These three columns signify “from”, “to” and “becomes” for our values. The first two, signal the range of values we want to classify. In this case, we don’t have wide ranges, we only want to reclassify one value (our NAs). But in order to do we still need to tell the classify() function what values, and their ranges, it is working with. Finally, the third column represents the value that the particular range of values we identified should become.

So, we can read that matrix as saying,“for cells with an NA value, convert them to 0s. Like wise for cells with a 0 value, convert them to 0s, and for cells with a value of 1, convert them to 1s.” Now, in this case, we aren’t converting our 0s and 1s to anything else, but that’s the basic logic. After we make our matrix, we can feed it into our classify() function after our rcl = argument.

elev_sub <- classify(elev_sub, rcl = elev_mat)

plot(elev_sub)

Our Reclassified Elevation Raster

So, what changed? Now, those blank or white cells are now 0s! This seems trivial, but, remember, NAs can sometimes wreck havoc on us when we are trying to run calculations, so converting to 0s will help us out further down in our workflow. This is not something you need to do every time you get NA in a raster, but you should consider it if you run into any weird calculation difficulties. Just makes sure that changing NA cells to a number does not fundamentally shift your analysis.

Pulling Forests and Parks

Next, we’ll bring in our landcover data and pull out the forests. We can do this using a simple == operator, but I want to show how we can reclassify our landcover dataset to get the same result. This is overkill, but I want you to see how can we can make use of this classification matrix approach in some different ways.

landcover <- rast(here::here("labs_exercises/Nov9_OR 2023-10-25_RasterIntro_DRAFT/Nov_9th In Class/Nov_9th In Class/Nov_9th/Raster_Intro/landcover/"))

forest_mat <- c(1, 0,
         2, 0,
         3, 1,
         4, 0,
         5, 0,
         NA, 0)

forest_mat <- matrix(forest_mat, ncol = 2, byrow = TRUE)

forests <- classify(x = landcover, rcl = forest_mat, others = 0)

Our Original Landcover Raster

Let’s break it down. We read in our landcover layer, and then we make our classification matrix. Now, note that this time we do not have a third “becomes” column, we just make use of the “from” and “to” columns. And we are reclassifying everything that is not forest land (category 3) to a 0 and classifying forests as 1. One additional bit here is that I also made use of the others = argument. This is a catch all argument that means that if my classification matrix is missing anything, to just convert that value to whatever we set for others =. In this case, if there is some unknown land classification we didn’t happen to catch we also convert those to 0.

Our Reclassified Forest Layer

Now, let’s work on our parks. First, we’ll read in our vector parks file using st_read() and then we will rasterize our vector layer by using the rasterize() function from terra. Note, this is different than what we did with our roads. For our roads, we kept it as a vector file, we just converted into a terra vector. But we’re going to want to work with our parks file with our other rasters which means it should also be a raster. This will be clear when we combine our layers. Just note that with rasterization we can offer a second raster for the parks object to share the same extent and resolution, and we convert all other values to 0s.

parks <- st_read(here::here("labs_exercises/Nov9_OR 2023-10-25_RasterIntro_DRAFT/Nov_9th In Class/Nov_9th In Class/Nov_9th/Raster_Intro/Parks.shp"))

Reading layer `Parks' from data source
  `/home/jamaalgreen/Documents/Courses/Fall_2024/cpln_503_fall_2024/labs_exercises/Nov9_OR 2023-10-25_RasterIntro_DRAFT/Nov_9th In Class/Nov_9th In Class/Nov_9th/Raster_Intro/Parks.shp'
  using driver `ESRI Shapefile'
Simple feature collection with 987 features and 12 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 2257090 ymin: 145086.3 xmax: 2475883 ymax: 360116.7
Projected CRS: NAD83 / Pennsylvania South (ftUS)

parks_rast <- rasterize(x = parks, y = elev, background = 0)

Our Nice Parks Raster

Bringing It All Together

Alright, we have the four criteria we need to get our final count of cells. Our final step is combining our layers and figuring out which cells meet all four criteria. To do that we can simply add our rasters together using map algebra!

rast_combined <- parks_rast + forests + roads_2mi + elev_sub
rast_combined <- mask(x = rast_combined, mask = elev)

plot(rast_combined)

Our Final Raster

Okay…we have one more new function here- mask(). We can use mask() to convert cells not covered by our mask layer to NAs. In this case, it’s helpful because we have many cells that are coded as 0s outside of our extent of interest and that will inflate the number of 0 cells. So, we convert those areas outside of the extent (in this case, the county) to NA.

Finally, how can we get a count of our cells for our categories? For that, terra gives us a handy function, freq(), which will give us the frequencies/counts of cells by some category. Let’s take a look.

freq(rast_combined)

  layer value   count
1     1     0  134733
2     1     1  748324
3     1     2 1520694
4     1     3  301379
5     1     4   34436

And with that, we’ve done a nice little bit of raster work. We have identified areas that meet our criteria based on both vector and raster data files, and we figured the number of cells that correspond to the number of conditions met. This is a less complicated example compared to other raster based workflows, but I hope it gives you a feel for the rhythm of raster analysis and highlights some of the ways that rasters may be useful for you in the future.

Again, if this is something that you are vibing with, then definitely think about “Modeling Geographic Space”.

Have an appropriately spooky Halloween